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Efficient  Numerical  and  Analog  Modeling 
of  Flicker  Noise  Processes 

J.   A.    Barnes  and  Stephen  Jarvis,    Jr. 

It  is  shown  that  by  cascading  a  few  simple  resistor  - 
capacitor  filters,    a  filter  can  be  constructed  which  generates 
from  a  white  noise  source  a  noise  signal  whose  spectral 
density  is  very  nearly  flicker,     |f  |      ,    over  several  decades 
of  frequency    f.     Using  difference  equations  modeling  this 
filter,    recursion  relations  are  obtained  which  permit  very 
efficient  digital  computer  generation  of  flicker  noise  time - 
series  over  a  similar  spectral  range.     These  analog  and 
digital  filters  may  also  be  viewed  as  efficient  approximations 
to  integrators  of  order  one -half. 

Keywords:    Analog  noise   simulation;   computer  noise   simulation; 
digital  filters;  flicker  noise;  fractional  integration;  recursive 
digital  filters. 

1.     INTRODUCTION 

One  of  the  most  pervasive  noise  processes  observed  in  electronic 
systems  is  flicker  noise,    whose  spectral  density  varies  as    |f  |         over 
many  decades  of  frequency,    often  beyond  the  low  frequency  limit  to 
which  the  spectral  density  can  reasonably  be  measured.      (See,    for 
example  [l]   -  [6]  and  their  bibliographies.  )    While  its  occurrence  is  very 
widespread,    its  cause,    or  causes,    is  as  yet  uncertain,    despite  the 
attention  of  many  investigators. 

In  order  to  study  flicker  noise  and  to  simulate  the  behavior  of 
systems  with  flicker  noise,    several  authors  have  presented  mathematical 
models  which  generate  a  discrete  [7],  [8]  or  continuous  [9]   -  [14]  flicker 
noise  from  a  white  noise  process.     Because  of  the  intrinsically  long 


correlation  time  associated  with  the  flicker  noise  process,    the  models 
have  required  large  computer  memories,    or  extensive  computation,    or 
both.     The  purpose  of  this  paper  is  to  present  an  algorithm  which  permits 
one  to  generate,    in  a  very  efficient  manner  with  negligible  computer 
memory  requirements,    a   sequence  of  numbers  with  a  very  nearly  flicker 
noise  spectral  density  over  several  decades  of  frequency.     It  is  derived 
from  a  cascade  of  simple  resistor -capacitor  filters  which,    when  realized, 
also  permits  one  to  generate  from  a  white  noise  source  an  analog  noise 
signal  which  has  a  spectral  density  which  is  very  nearly  flicker  over 
several  decades  of  frequency. 

2.      THE  FILTER  CASCADE 
Consider  the  filter  shown  in  figure  1.     Its  voltage  transfer  function 

is: 

g(C0)  =  YM=      i+jq>T8 

s^'      x(co)      1  +  jwOi  +  xa)  [L> 

where  Tx  =  R:C  (2) 

t3=R2C  (3) 

and  C0=  2  77f.  (4) 

Equation  (1)  is  valid  provided  the  loading  on  the  output,    y,    of  the  filter 
is  negligible.     This  assumption,    thus,    might  require  the  filter  to  be 
followed  by  an  isolation  amplifier  in  practice.     As  it  will  be  seen, 
however,    it  is  possible  to  construct  a  complete  flicker  filter  with  only 
one  amplifier  and  still  not  compromise  the  validity  of  eq  (1). 

For  low  frequencies   (<jO(Ti  +  t2)  <<  l ),    one  finds  that  g(C0)  «  1;  and 

T2 

for  high  frequencies   (to  t2  >>1),    one  finds  g(C0)  « ; .     If  one  were  to 

T*i  +  Ts 

have  a  sequence  of  such  filters  such  that  for  the  i-th  filter 

A{)  =  (Pf^K     i=  1.2....N,  (5) 

where    jS   is  some  constant  such  that  (8  <  1, 
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T2(i)   =    (/3)1  T  8<°> 


(6) 


and  if  these  filters  were  cascaded  with  appropriate  isolation  amplifiers 
between  successive  stages   (to  eliminate  loading  effects),    the  overall 
voltage  transfer  function  would  become: 


N[     1+   iC0T2(°V  1 

G(co)=n    — — ■""  a    p         .    . 

i=lLl  +jco(tx(o)  +  rJo))^ 


(7) 


One  may  consider  a  set  of    n    sections  of  this  filter  cascade 
corresponding  to    i  =  j  +  1,    j  +  2,    .  •  .  j  4-  n,    where    j    is  arbitrary. 
This   subset  of  filters  would  change  the  (magnitude)  transfer  function  by 
the  factor 

/  r  (o)      \n 
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in  the  frequency  interval 
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If  this  filter  is  to  approximate  a  filter  whose  transfer  function  varies 


as 


'G*(CO)|2=A8|cor 


where  A  is  a  constant,    then 
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-  1  in  eq  (8)  and,    thus,    eq  (11)  becomes 

A  more  precise  approximation  to  a  flicker  filter  is  obtained  by  increasing 

/8  (i.e.,    the  ratio  of  the  t's  given  in  eq  (12)).     However,    the  frequency 

N 
range  of  the  complete  filter  is  proportional  to  P       .      Thus,    to  improve 

the  precision  without  reducing  the  frequency  range  both  $   and  the 

number  of  sections,    N,    must  be  increased. 

It  is  of  value  to  take  as  a  practical  example  the  particular  values: 

N  =  4  (13) 

Ti(o)  +  Ta(o)=  3Ts(o)  (u) 


and  use  the  variable 


jO)T2(l).  (15) 


The  filter's  transfer  function  then  becomes  (from  eq  (7)) 

1 

—  p 

This  transfer  function  should  be  expected  to  have  approximately  an  s 

behavior  over  a  relative  frequency  range  of    (8       =  9   ,    nearly  four 

decades  of  frequency.     In  fact,    the  magnitude  error 

i 
e=   |G(W)|    -  A|  s"2|,  (17) 

l 
is  within  ±  -g  dB  relative  to  A  |  s    '        over  the  range  0.  2  5  ^    |s  |   ^  1000,    a 

relative  frequency  range  of  4000  where  A  =  0.  43.     The  location  of  the 

approximated  frequency  band  depends  on    T       .     For  the  value  T       =  0.  04 

2  3 

seconds,    the  approximating  band  is   (1  Hz,    4  kHz);  the  error  is  shown  in 
figure  2. 

One  may  forget,    for  the  time  being,    the  heuristic  arguments  which 
led  to  generating  eq  (16)  as  an  approximation  to  the  function 


1 
G*(CO)  =  As"2.  (18) 

Equation  (16)  simply  gives  a  mathematical  function  which  approximates 
eq  (18)  rather  well  over  a  substantial  frequency  range.      Thus,    one  could 
consider  an  impedance,    zi  (CO),    defined  by 

zi(to)  =  RoG(CO),  (19) 

where  the  constant    Ro     has  the  dimensions  of  ohms.      The  exact  realiza- 
tion of    zi(C0)    may  be  obtained  by  a  partial  fractions  expansion  of  eq  (16) 
and  one  obtains  the  impedance  shown  in  figure   3   for  R0   =  20kO      and 
T21  =  0.  04  seconds.      (The  exact  realization  of  an  impedance  given  by 
eq  (18)  would  be  much  more  difficult.  )      Thus,    figure   3  is  an  approxi- 
mation to  a  "fractional  capacitor"  [  15]   -  [17]  of  order  one -half. 

This  impedance  can  be  used  in  a  straightforward  fashion  in  con- 
nection with  a  single  operational  amplifier  to  realize  a  transfer  function 
given  by  eq  (16).     It  is  no  longer  necessary  to  consider  the  filter  cascade 
with  its  isolation  amplifiers  mentioned  above.     This  new,    active  filter 
is  realized  as  follows: 

An  operational  amplifier  is  a  device  whose  transfer  function  is 
approximated  by  a  large,    real,    negative  number,      -k,    over  a  large 
frequency  range  which  is  assumed  to  include  the  frequency  domain  of 
interest.     If  one  connects  an  operational  amplifier  with  two  impedances 
zx(co)    and  Z2(co)    as  shown  in  figure  4,    the  overall  transfer  function 
becomes 
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Thus,    if  one  chooses     z2    to  be  a  pure  resistance  of  value    Ro,    the  over- 
all transfer  function  is  just  G(Cd)    given  by  eq  (16).     Such  a  device  has, 
in  fact,    been  built  and  tested  at  the  National  Bureau  of  Standards. 

It  is  interesting  to  note  that  if  one  interchanges  the  impedances 
zi    and    Z2    in  figure  4,    one  obtains  a  prewhitening  filter  for  flicker  noise 
in  the  frequency  range  from  1Hz  to  4  kHz. 

3.     EFFICIENT  GENERATION  OF  APPROXIMATE 
FLICKER  NOISE  NUMBERS 
From  the  impulse  response  function 

Mt)  =  ^=  (22) 

corresponding  to  the  transfer  function 

one  can  generate  a  flicker  signal    eo(t)    from  a  white  noise  source    e.(t) 
by  the  convolution  integral  [18]: 

e,(T) 


e0 


(t)=r.dT^T-  (23) 


Discrete  flicker  numbers  can  be  approximated  by  a  discrete  analog 

convolution,    setting    t     =  n  6t: 

n 

m 

(eo)       =y    h(t     )(e.)  ,  (24) 

m      Z_,         m      i  m-n 

n=l 

where  the  (e.)       are  uncorrelated  random  deviates.     To  obtain  a  flicker 
in 

spectral  density  over  a  relative  frequency  range     r,    the  memory  length 
M   must  be  very  large   and  so  also  the  number  of  remembered  deviates 
(e.)       for  each  output     (eo)    .      This  is  not  a  very  efficient  use  of 
computer  memory  space  nor  of  computer  time. 

A  very  efficient  flicker  number  generator  can  be  constructed  from 
the  filter  cascade  introduced  in  section  2.     The  differential  equation  for 


10 


the  network  of  figure  1  is: 

y(t)  +  y(t)(Ri  +  R2)C  =  R2Cx(t)  +  x(t).  (25) 

This  can  be  approximated  by  the  difference  equation  for  the  discrete 


values    y     =  y(n  fit): 
n 


=   (1  6t  )  |       Ts 

Yn+1        [         Ta  +  T2  ;  yn       TX  +  t2  ~n+l  '    Ta   +  t2  ~n' 


,    6t   -  T5 
x    .  .  +  ■ x 


(26) 


with,    as  before,      T.  =  R.C.     Let 


and 


y  = 


R 


fit 


Ti  +   T2     ' 


Tl    +   T3 
The  recursion  relation  for    y       is  then 


(27) 


(28) 


yn+l  =  (1   "  yK  +  Rxn+1   "  (R  -  yK> 


(29) 


for  an  unloaded  filter.     As  in   section  2,    we  shall  again  consider  four 
such  filter  sections,    with  parameters    y       and  R      ,    i  =   1,    2,    3,    4,    inputs 
x  and  outputs    y        .     Cascading  the  filters,    the  output  of  filter    i    is 

the  input  of  filter    i  +  1: 


y  «   =  x  (2> 

n  n 


(2) 


(3) 


y  ^   -  x  (4) 

n  n 


(30) 


Note  that  this  is  equivalent  to  assuming  the  existence  of  isolation  ampli- 
fiers between  sections  of  the  analog  filter.      That  is,    eq  (2  9)  gives  the 
exact  filter  section  response  (no  loading  effects)  and  its  output  is  exactly 
the  input  to  the  next  stage   (eq  (30)). 
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Then,    'with  the  earlier  definitions: 


(i) 


R(l)  =  R  =  1/3  , 


and 


0   = 


(i) 


(i)   ,    ,(i) 


=   (1/3)* 


'1        +    T"2 

In  this   situation  the  higher -numbered  filters  have  the  shorter  time 
constants.     It  seems  reasonable  to  let  the  shortest  time  constant,    Ta' 
be  equal  to  the  inverse  of  the  Nyquist  frequency,    f 


NQ       (2)(6t) 


Then 
and 


fit  =   l-  ^ 


,(i)  =  I(I\9"21    (i=  1.4)    • 


(31) 
(32) 
(33) 

(4) 
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(34) 
(35) 


The  resulting  recursion  relations  for  the  four -filter  cascade  become: 
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(3)  (3)_    /485\    ■  W+Ix  0) 

Xn+i       yn+i~    \486/    Yn     +  3     n+i 


1457  \      (1) 
4374/   Xn 

161  \       (3) 


186 


/       n 


CW    =y   (3>=     /«)     y(3)+Ix(3)    -     fll)       x(l 
n+i       Yn+i        I  54/     Yn     +  3     n+i         154/      xn 


(s) 


(4)  ,1        (4)        1       (4) 

3     n+i       b     n 


4)       /5\      (., 


(36) 


When  the    x        are  uncorrelated  random  deviates,    and  the  Ix      ,    y       J 
n  I    o  0     I 

are  taken  to  be  zero  initially  (to  minimize  the  transient),    the  above  recur- 
sion equations  yield  a  set  of  numbers,    y  4,    which  approximate  flicker 
noise  over  a  range  of  about  1000  periods  of  6t.    (The  addition  of  a  fifth, 
lower  frequency,    or  i  =  0,    stage  should  extend  this  range  to  nearly  10,  000 
periods.)     The  recursion  is  neatly  accomplished  on  a  digital  computer 
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with  a  random  number  generator  and  storage  locations  for  the  eight 

,.    ,.  „.    .  ,        .,  ,  (s)       (3)        (4)         (1)         (2) 

distinct  coefficients  plus  the  seven  numbers    x     ,    x      ,    x       ,    y       ,    y        , 

r  n         n  n  n  n 

y (3),  y  U>    - 

n  n 

Figure  5  shows  a  Fortran  program  written  for  the  generation  of  1024 

numbers  which  approximate  a  flicker  noise  sample  using  eqs   (36). 

Figure  6  shows  a  sample  of  N  =   1024  numbers  obtained  using  this  program. 

This  block  of  data  was  subjected  to  a  time -domain  statistical  analysis 

developed  by  Allan  [19],    [20].     Figure  7  shows  the  Allan  variance,    (7(t), 

for  independent  samples.     For  pure  flicker  noise  of  arbitrary  length, 

a(T)  has  no  T -dependence.     The  deviations  of  the  plot  from  a  straight 

horizontal  are  probably  due  to  the  small  sample  size. 

Figure  8  is  a  plot  of  the  actual  impulse  response  function  of  the 
digital  filter  described  in  the  Fortran  program  of  figure  5. 

4.     ACKNOWLEDGMENT 
The  authors  wish  to  thank  Mr.    D.   W.  Allan  of  the  National  Bureau 
of  Standards  for  his  assistance  in  the  statistical  analysis  of  the  data. 


13 


PROGRAM  FLICKER 
DIMENSION  V(5,  2) 
C  SET  INITIAL  VALUES  V(I,  1)  TO  ZERO  MEAN  VALUES 

DO  1     1=1,  5 

1  V(I,1)  =  0. 

C  GENERATE  AND  PRINT  N  FLICKER  NUMBERS 

N  =   1024 
DO  2  1=1,  N 

C  SELECT  RANDOM  V(l,  2)  UNIFORMLY  DISTRIBUTED  ON  (-|,  |) 

V(l,  2)  =  RANF   (-1)   -0.  5 

C  SOLVE  RECURSION  RELATIONS 

V(2,  2)  =  .  999771*V(2,  1)  +  .  333333*V(1,  2)  -  .  333105*V(1,  1) 
V(3,  2)  =  .  997942*V(3,  1)  +  .  333333*V(2,  2)  -  .  331276*V(2,  1) 
V(4,  2)  =  .  981481  *V(4,  1)  +  .  333333*V(3,  2)  -  .  314815*V(3,  1) 
V(5,  2)  =  .  833333*V(5,  1)  +  .  333333*V(4,  2)  -  .  166667*V(4,  1) 
PRINT  11,  V(5,  2) 
11      FORMAT  (2XF12.6) 

C  RESET  V 

DO  3  J=l,  5 
3     V(J,  1)=V(J,2) 

2  CONTINUE 
CALL  EXIT 
END 


Figure  5.     Fortran  program  for  computing  N  approximate 
flicker  numbers  by  recursion  relations,    eq  (36). 
("N  =  1024"  has  been  specified.  ) 
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